#include <mkl.h>
#include "CPPDiffusion2DADI.h"
      
double CPPDiffusion2DADI::BoundaryCondition() {
	return m_boundaryCondition;
}

void CPPDiffusion2DADI::BoundaryCondition(double boundaryCondition) {
	m_boundaryCondition = boundaryCondition;
}

int CPPDiffusion2DADI::NumberOfNodesInX() {
	return m_nx;
}

void CPPDiffusion2DADI::NumberOfNodesInX(int nx) {
	m_nx = nx;
}

int CPPDiffusion2DADI::NumberOfNodesInY() {
	return m_ny;
}

void CPPDiffusion2DADI::NumberOfNodesInY(int ny) {
	m_ny = ny;
}

double CPPDiffusion2DADI::DeltaX() {
	return m_deltax;
}

void CPPDiffusion2DADI::DeltaX(double deltax) {
	m_deltax = deltax;
}

double CPPDiffusion2DADI::DeltaY() {
	return m_deltay;
}

void CPPDiffusion2DADI::DeltaY(double deltay) {
	m_deltay = deltay;
}

double CPPDiffusion2DADI::DeltaT() {
	return m_deltat;
}

void CPPDiffusion2DADI::DeltaT(double deltat) {
	m_deltat = deltat;
}

double CPPDiffusion2DADI::Alpha() {
	return m_alpha;
}

void CPPDiffusion2DADI::Alpha(double alpha) {
	m_alpha = alpha;
}

void CPPDiffusion2DADI::NextStep(double * u) {

	  double * u_half = new double[m_nx * m_ny];

      for (int j = 1; j < m_ny - 1; j++) {
        StepI(u, u_half, j);
      }

      for (int i = 1; i < m_nx - 1; i++) {
        StepII(u, u_half, i);
      }

	  delete[] u_half;

}

void CPPDiffusion2DADI::StepI(double * u, double * u_half, int j) {

	int n = m_nx;
	int nrhs = 1;
	double * d = new double[m_nx];
	double * dl = new double[m_nx - 1];
	double * du = new double[m_nx - 1];
	double * b = new double[m_nx];
	int ldb = m_nx;
	int info = 0;

	d[0] = 1.0;
	du[0] = 0.0;
	b[0] = m_boundaryCondition;

	double rx;
	double ry;

	for (int i = 1; i < m_nx - 1; i++) {
		rx = (m_alpha * m_deltat) / (m_deltax * m_deltax);
    ry = (m_alpha * m_deltat) / (m_deltay * m_deltay);
		
		dl[i - 1] = - rx;
		d[i] = 1 + 2 * rx;
		du[i] = - rx;
		b[i] = ry * u[(j - 1) * m_nx + i]
		     + (1 - 2 * ry) * u[j * m_nx + i]
			   + ry * u[(j + 1) * m_nx + i];
				
	}

	d[m_nx - 1] = 1.0;
	dl[m_nx - 2] = 0.0;
	b[m_nx - 1] = m_boundaryCondition;

	dgtsv(&n, &nrhs, dl, d, du, b, &ldb, &info);

	for (int i = 0; i < m_nx; i++) {
		u_half[j * m_nx + i] = b[i];
	}

	delete[] d;
	delete[] du;
	delete[] dl;
	delete[] b;
}

void CPPDiffusion2DADI::StepII(double * u, double * u_half, int i) {

	int n = m_ny;
	int nrhs = 1;
	double * d = new double[m_ny];
	double * dl = new double[m_ny - 1];
	double * du = new double[m_ny - 1];
	double * b = new double[m_ny];
	int ldb = m_ny;
	int info = 0;

	d[0] = 1.0;
	du[0] = 0.0;
	b[0] = m_boundaryCondition;

	double rx;
	double ry;

	for (int j = 1; j < m_ny - 1; j++) {
		rx = (m_alpha * m_deltat) / (m_deltax * m_deltax);
    ry = (m_alpha * m_deltat) / (m_deltay * m_deltay);
		
		dl[j - 1] = - ry;
		d[j] = 1 + 2 * ry;
		du[j] = - ry;
		b[j] = rx * u_half[j * m_nx + (i - 1)]
		     + (1 - 2 * rx) * u_half[j * m_nx + i]
			   + rx * u_half[j * m_nx + (i + 1)];
				
	}

	d[m_ny - 1] = 1.0;
	dl[m_ny - 2] = 0.0;
	b[m_ny - 1] = m_boundaryCondition;

	dgtsv(&n, &nrhs, dl, d, du, b, &ldb, &info);

	for (int j = 0; j < m_ny; j++) {
		u[j * m_nx + i] = b[j];
	}

	delete[] d;
	delete[] du;
	delete[] dl;
	delete[] b;

}
